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Fixed-Node Diffusion Monte Carlo of Lithium Systems 


K. M. HasclQ and L. MitatQ 
Physics Dept., North Carolina State University 

We study lithium systems over a range of number of atoms, e.g., atomic anion, dimer, metallic 
cluster, and body-centered cubic crystal by the diffusion Monte Carlo method. The calculations 
include both core and valence electrons in order to avoid any possible impact by pseudopotentials. 
The focus of the study is the fixed-node errors, and for that purpose we test several orbital sets 
in order to provide the most accurate nodal hypersurfaces. We compare our results to other high 
accuracy calculations wherever available and to experimental results so as to quantify the the fixed- 
node errors. The results for these Li systems show that fixed-node quantum Monte Carlo achieves 
remarkably high accuracy total energies and recovers 97-99 % of the correlation energy. 

PACS numbers: 02.70.Ss, 31.15.V-, 71.10.-w, 71.15.-m, 71.15.Nc, 71.20.Dg 
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I. INTRODUCTION 

Quantum Monte Carlo (QMC) methods have been ap¬ 
plied to a great variety of electronic structure problems 
over the past three decades. These calculations have pro¬ 
vided a number of high accuracy results for properties 
such as cohesion and binding energies, excitations, reac¬ 
tion barrier heights, defect formation energies, and other 
quantities; they are t ypical ly in excellent agreement with 
available experiments®^. In addition, the calculations 
have shed new light on correlation effects in various sys¬ 
tems, and therefore have become valuable as benchmarks 
for other methods and comparisons. The most impor¬ 
tant strength of this approach is that the many-body 
Hamiltonian is employed directly, and thus the electron- 
electron interaction and particle correlations are treated 
explicitly in a many-body manner. Additionally, QMC 
methods can be applied to large systems so that prop¬ 
erties of solids can be calculated by using supercells and 
extrapolations to the thermodynamical limit. 

Diffusion Monte Carlo (DMC) projects out the ground 
state of a system by applying the projection operator 
exp(— tH), where H is the Hamiltonian, to the trial wave 
function 'L'p- In the large imaginary time limit r —> oo, 


the ground state of a given symmetry is obtained. One 
of the key fundamental limitations to achieving exact 
results is the so-called fixed-node approximation which 
en ables one to avoid the well-known fermion sign prob- 
len J 19 l 39 l . The fixed-node approximation is very difficult 
to improve upon since the corresponding energies are typ¬ 
ically very small, e.g., a few percent of the correlation en¬ 
ergy where the correlation energy is itself a small fraction 
of the total energy. Therefore systematic improvements 
of the nodes through minim ization of the total energy or 
variance of the energjPUM f or a given trial wave func¬ 
tion is laborious and often very costly. Alongside the ob¬ 
viously desirable achievements needed in improving the 
algorithmic efficiency, in increasing the speed and qual¬ 
ity of wave function optimizations, and in exploring new 
functional forms for orbitals and wave functions, a signif¬ 
icant amount of work remains to be done in simply un¬ 
derstanding the root origin of the nodal errors by system¬ 
atically quantifying the dependence of energy biases on 
the nodal defects. Recently, we have analyzed the impact 
of the electron density on the nodal bias in a set of free 
atom/ionic systems, and we found the fixed-node errors 
being proportional to the density in this particular class 
of s ystem s for both spin-polarized and spin-unpolarized 
cased®®. More insights have been gained on the role of 
the basis setd® i n the error although the findings were 
not always easy to utilize in different system^!. 

As a testbed to expand this study for more compli¬ 
cated cases, a series of lithium systems is attractive for 
several reasons. First, inclusion of the core electrons in Li 
calculations is computationally feasible. That enables us 
to avoid any additional, more complicated analysis that 
is necessary whenever pseudopotentials or effective core 
potentials are being employed. This ensures that any 
missing amount of the binding energy, cohesive energy, 
or correlation energy is caused solely by the fixed-node 
approximation. Second, the exact wave function for a 
free Li atom has a relatively simple nodal surface which 
is already well approximated at the Hartree-Fock level. 
In addition, small Li systems have been studied with the 
FN-DMC method before. In this work we significantly 
expand upon previous studies with calculations of Li 4 
cluster and Li solid in its equilibrium body-centered cubic 
structure. By comparing our fixed-node diffusion Monte 
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Carlo (FN-DMC) results with other accurate calculations 
and experimental results corrected for zero-point motion 
(DMC is carried out in the Born-Oppenheimer approxi¬ 
mation), we can assess the magnitude of the fixed-node 
errors with high accuracy. These systems represent a 
variety of environments for the bonding and include di¬ 
rectional bonds, multi-center bonding, and delocalized 
metallic bonding. It is therefore an interesting question 
to understand how the fixed-node bias changes once Li 
enters bonding in the setting of molecular bonds or pe¬ 
riodic boundary conditions. Based on these results we 
establish a systematic picture of the nodal errors in Li 
systems and corresponding accuracies for energy differ¬ 
ences. 


II. METHODS 

In DMC, we solve for the ground state solution of 
Schrodinger’s equation 

T 0 = lim exp{—rJf}$T (1) 

T—>■ OO 

where is the Born-Oppenheimer Hamiltonian. The 
antisymmetric nature of fermion systems poses a chal¬ 
lenge to the naive application of the DMC algorithm and 
leads to the fermion sign problem 3 39 52 . This is because 
for a given boundary value problem, the eigenstate with 
the lowest eigenvalue will be a symmetric state. In light 
of this well-known difficulty, perhaps the simplest and 
most straightforward way to circumvent the sign prob¬ 
lem is the fixed-node approximation. 

Under the fixed-node approximation, we impose a 
boundary condition at the nodes of the trial wave func¬ 
tion and maintain them for the duration of the simula¬ 
tion. The nodes form a hypersurface defined implicitly 

by 

0 = {R;T t (R) =0}. (2) 

The assumed nodal hypersurface creates boundaries that 
constrain the solution in each nodal cell and preserve the 
overall fermionic anti-symmetry of the total wave func¬ 
tion, preventing thus any appearance of “signs.” This 
allows one to ignore the sign of the wave function inside 
the nodal cell and to carry out the DMC algorithm within 
each nodal cell 

TtTo > 0. (3) 

Unfortunately, doing this exactly is a tall order as it re¬ 
quires that for an N electron system one must have a 
description of the exact (31V — l)-dimensional lrypersur- 
face Q. Solving for such a hypersurface directly is beyond 
our means, and instead we proceed by using nodal sur¬ 
faces from approximate wave functions. Because we use 
a nodal hypersurface that is not exact, the solution will 
have a higher energy than the exact ground state, i.e., 
the total energy computed via FN-DMC is a variational 
upper bound to the exact energy 3 ®. Further details on 
the FN-DMC method can be found elsewhere^. 


The trial functions used in this study are of the Slater- 
Jastrow type 

*t(R) = XI Cfedetlbildettl^] exp(U) (4) 

k 

where the one-particle orbitals are obtained from 
Hartree-Fock (HF) or Density Functional Theory (DFT). 
More details on the trial functions and their optimiza¬ 
tions, Jastrow correlation factors and DMC calculations 
for periodic boundary conditions can be found in the re¬ 
cent review 27 . 


III. RESULTS 

III.l. Lithium atom & electron affinity 

In order to calculate the total energy of the 
lithium atom, we use a restricted Hartree-Fock, single¬ 
determinant wave function. Since at the HF level the 
spin-up and spin-down subspaces are independent, the 
minority spin channel contains only one electron, and the 
wave function’s node exists in only in the spin majority 
subspace. We can visualize this subspace of the nodal hy¬ 
persurface by considering the wave function when 1 and 2 
label the electrons in the same spin channel. Then it fol¬ 
lows from the form of the HF determinant that the node 
is given by the condition rq = r 2 - The electron 1 there¬ 
fore “sees” the node as a sphere which passes through 
the position of electron 2 and is centered around the nu¬ 
cleus. The wave function will be equal to zero if electron 
1 occupies any point on the spherical nodal surface. 

The exact HF nodal hypersurface in the full 6D space 
is a 5D hyperboloid given by the implicit equation 
x i + y\ + Z? = x 2 + V 2 + z 2 - As pointed out by Still- 
inger et alW- 1 , this is not strictly exact, as the correlation 
with the electron in the spin-down channel will cause de¬ 
formations away from a perfect sphere. For example, the 
excitation 2s 1 2p 2 will have a contribution to the exact 
ground state and would in principle lead to a departure 
from the single particle node (i.e., the sphere will slightly 
deform to ellipsoid or perhaps a more complicated sur¬ 
face that would depend on the position of the minority 
spin electron). It is therefore quite remarkable that the 
HF nodal surface seems to be so accurate: the total en¬ 
ergy with the HF nodes, shown in Table [TJ is accurate 
to « 0.05(1) mHa and the fixed-node bias is less than 
0.1 % of the correlation energjffD. This points out that 
the correlation is basically completely captured by the 
Jastrow-like effect and it affects the 5D hyperboloid only 
marginally. (This contrasts with the Be atoms where the 
nodal surface is strongly affected by correlation:!^!) Our 
own calculated fixed-node error in the single atom energy 
is ss 0.05(1) mHa, much smaller than, say, the chemical 
accuracy (« 1.6 mHa). This also suggests that any fixed- 
node errors in the aggregate species from Li atoms will 
be essentially identical to the error in the binding or co¬ 
hesive energy. 
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TABLE I. Comparison of theoretical results for the total en¬ 
ergy of a lithium atom where FN-DMC energy has been ob¬ 
tained with the HF nodes. 


source total energy (Ha) 

FN-DMC present work -7.47801(1) 
FN-DMC Bressanini et a/P® -7.47803(5) 
exaclP -7.47806032 


TABLE II. Comparison of the latest calculation and measure¬ 
ment with FN-DMC results for the electron affinity for Li (in 
Ha). 


Author Method EA (Ha) 

FfschepSI extrap. MCHF 0.022698 

present work FN-DMC HF single det. 0.0201(1) 
present work FN-DMC 2 configs. 0.02279(5) 

Haeffled^ expt. _0.0227129(8) 


Using the value for the ground state total energy of the 
Li - ion from reference®!, which has a fixed-node error 
comparable to the neutral atom, we can compute the 
electron affinity of a Li atom. The electron affinity is 
given by 

EA(Li) = -Eo(Li) — E 0 (Li _ ). (5) 

Lithium has a positive electron affinity, meaning the an¬ 
ion is more stable than the neutral atom. The HF limit 
of the total energy of Li” has been computed® to be 
— 7.4282320 Ha so that in the HF approximation the ad¬ 
ditional electron is not bound. As it is well-known, the 
correlation effects are crucial for describing the electron 
affinity with accuracy comparable to experiment I 24 -l 

Using a wave function composed of two configuration 
state functions, i.e., HF plus 2s 2 —1- 2 p 2 excitation (a 
symmetry adapted linear combination of determinants), 
given the data from our recent calculation^®! for the 4 
electron ion yields an FN-DMC electron affinity with an 
excellent accuracy compared to experimental measure¬ 
ment 21 . A summary of the best theoretical and exper¬ 
imental data are compared with this work in Table [TT] 
In terms of the fixed-node error, while the two species in 
the calculation share geometric details (central potential 
in free space), the nodal hypersurface changes when the 
4th electron is added and the anion shows similar nodal 
shape as the isoelectronic Be atom. The poor quality of 
the result for the electron affinity using only a single de¬ 
terminant trial function for Li” stands to illustrate that 
the extreme accuracy of the RHF nodal hypersurface for 
3 electrons is not typical and rather results from a fortu¬ 
itous coincidence. 


Lii 



FIG. 1. Schematic depiction of the D 2 h Li 4 parameters 


III.2. Li 2 

The Li 2 dimer is a more complicated system, with 
increased number of electrons while also changing the 
overall real-space geometry from one central potential 
with spherical symmetry to a cylindrical symmetry. The 
nodes of the lithium dimer have been a studied num¬ 
ber of times 12 17 . The best result in the literature (in 
reference 12 -!) exceeds the quality of our single configura¬ 
tion result, -Efn-dmc = —14.9932(7), which has almost 
2 mHa fixed-node error. The fixed-node error of the best 
wave function in the literature is « 0.2 — 0.3 mHa, 
with a total energy of —14.9952(1). Using the value for 
the HF energy reported by Filippi and Umrigad 17 !, this 
recovers « 99.8% of the correlation energy. Bressanini 
et aZ! l2j 'have pointed out that a five configuration wave 
function basically gets the nodal surface close to exact, 
and one can get accurate energies at the fixed-node ap¬ 
proximation level with errors of the order of the order 
of 0.1 % as well. Interestingly, only the three lowest ex¬ 
citations were really involved, 2er 2 —► l7r 2 , 2cr 2 —► 2cr 2 , 
2ag —» 3cr 2 . Since this system was studied exhaustively, 
we did not repeat the muti-reference calculation, and in¬ 
stead we quote the results of Bressanini et a/P® hereafter. 


III.3. Li 4 

Li 4 and its properties have been studied by sev¬ 
eral metho ds includi ng the basis set correlated ap¬ 
proached 9 23 32 37 l 4 ° l 47 l The most stable configuration 
of four Li atoms is a molecule with D 2 h symmetry, a 
planar rhombus geometry, and singlet electronic ground 
stat e! 9 ! 23 ! 37 ! 40 ! 47 !. The geometry of the D 2 h Li 4 is depicted 
schematically in Figure [l] This can be understood as the 
result of a Jahn-Teller distortion of the more symmet¬ 
ric geometry of a square!®]. Li 4 exhibits a “three-center” 
bonding pattern where two electrons are shared inside 
each of the two triangles formed by bisecting Figure [l] 
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TABLE III. Fixed-Node DMC total energy for trial wave func¬ 
tions from different levels of theory testing unoptimized nodal 
surfaces for use as DMC trial wave-functions. 


theory 

size of virtual space 

(Ha) 

RHF 

0 

-30.0177(5) 

CI-SD 

9 

-30.0184(5) 

CI-SDTQ 

9 

-30.0174(4) 

CI-SD 

15 

-30.0228(4) 

CI-SDTQ 

15 

-30.0179(4) 

CI-SD 

19 

-30.0162(4) 

CI-SDTQ 

19 

-30.0179(6) 


TABLE IV. FN-DMC results for different basis sets with trial 
wave functions from CI-SD calculations using 15 virtual or¬ 
bitals and then optimized in with respect to VMC total en¬ 
ergy. 


basis set total energy (Ha) 

Roos Aug. DZ ANO (4s3p2d) -30.02127(5) 

Roos Aug. TZ ANO (4s4p3d2f) -30.02119(6) 
aug-cc-pCVTZ (7s6p4d2f) -30.02263(6) 


along the “r^” lind^l 

In Table m we compare the nodes of the SCF wave 
functions (unoptimized) for different levels of Cl in order 
to illustrate the behavior of such expansions and to select 
the best starting place for our QMC trial wave function. 
For each level of Cl, we used complete expansions but 
limited the number of virtual orbitals in the active space. 
It is clear that the nodes do not improve systematically 
for larger active space and higher level of theory as the 
Cl total energies do. 

Since FN-DMC errors associated with the basis set are 
not very systematic, we also tested the nodal surfaces of 
several basis sets to minimize these errors. Although not 
fully complete, the results seem to support the conclusion 
of Bressanini et al. that for Li systems saturating the s 
channel is more important than adding additional high 
angular momentum basis functions. These results are 
listed in Table |IV| After some initial testing of basis and 
multi-determinant expansions, we employed wave func¬ 
tions constructed from the aug-cc-pCVTZ basis and in¬ 
cluded the 15 lowest lying virtual orbitals into the CI-SD 
calculation. We re-optimized the weights of the result¬ 
ing 93 configuration state functions in the Cl expansion 
with VMC total energy minimization using a Levenberg- 
Marquardt algorithm. The geometry parameters have 
been computed a number of times in the literature, as 
reported in Table [V] organized by the value for the Li-Li 
distance labeled rc in Figure[lJ The results indicate that 
there is a short Li-Li bond ss 2.64 A (2r^ in Figure Jl]) 
and a longer Li-Li bond ss 2.99 — 3.0 A (rc in Figure IlF. 

The experimental and theoretical binding energies of 
Li 4 are given in Table VI Because the presence of the Jas- 
trow factor will influence the optimization of the multi¬ 


Time-Step Error for VMC optimized Wave Functions 



FIG. 2. The DMC energy extrapolated to t = 0 for a Lb 
molecule. Both wave functions originate from the same Con¬ 
figuration Interaction CI-SD calculation, but have had the de¬ 
terminant weights and Jastrow coefficients re-optimized with 
different assumed cut-off distances in the Jastrow, as de¬ 
scribed in the text. 


determinant expansion, it is not clear a priori what form 
of the Jastrow factor is optimal. In particular, since the 
bonds are not very strong and the bond lengths are some¬ 
what larger than in typical single-bonded situations we 
tested the range of the Jastrow cutoff distance parameter. 
We optimized the Jastrow coefficients and determinant 
weights for two different electron-ion Jastrow distances. 
The wave function with the so-called “short-range” Jas¬ 
trow effects had its electron-ion and electron-electron-ion 
terms extend to 2.45 bohr from each atom, i.e., to just 
less than half the smallest Li-Li distance. For the “long 
range” Jastrow, the electron-ion terms were allowed to 
extend to 7.5 bohr. The qualitative difference between 
these two Jastrows is that terms from different atoms in 
the short-range Jastrow do not overlap in the region oc¬ 
cupied by the three-center bonds. This difference in de¬ 
scription of the wave function translates into a difference 
in computational effort. In the “short-range” case, each 
electron will have non-zero three-body Jastrow terms as¬ 
sociated with only 1 atomic center at any given time; 
whereas in the “long-range” case, electrons have non-zero 
contributions to the Jastrow coming from each of the Li 
atoms surrounding the three-center bonding region. The 
resulting effect in the total energy of the wave function 
was only « 0.00013 per atom. 

We carried out a time-step extrapolation for both wave 
functions, shown in Figure[2j to ensure that the time-step 
error is <~ 0.05 mHa in the total energy, or an order 
of magnitude smaller than the statistical error bars in 
the binding energy. After correcting for the zer o p oint 
motion, which was about 3.12 mHartree per atonPZl, we 
find a binding energy of 0.723(3) eV. Note the reason¬ 
ably good agreement between the results although the 
basis set correlation methods did not include any corre¬ 
lation of the core (Is) states. This points out that the 
core states are already quite deep and do not affect the 
nodal surfaces significantly. One of the reasons is that 
any excitation which would correlate the Is level would 
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TABLE V. Summary of the optimized geometry parameters of D 2/1 Lb tested in this work, and the FN-DMC total energy for 
each. The trial wave function is an VMC energy optimized CI-SD expansion with 93 CSFs. 


author method ta (A) tb (A) rc (A) (Ha) 

RajP 3 DFT 1.298 2.759 3.050 -30.02132(2) 

Rousseau and ManP QCISD / CCSD(T) 1.323 2.700 3.007 -30.02155(2) 

Verdicchio et atP CCSD(T) 1.323 2.680 2.989 -30.02158(2) 

Wheeler et aZP CCSD(T)1.316 2.681 2.987 -30.02154(2) 


TABLE VI. Binding energies of LU, uncorrected for zero- 
point motion, are given in units of eV per atom. 


TABLE VII. Results of FN-DMC calculations for the F-point 
wave function of an 8 atom supercell comparing the nodal 
quality of orbital sets generated by DFT functionals. 


author 

method 

binding 



Alikhani and ShailP 

DFT 

0.61 

Exchange 

Correlation E to t (Hartrees) 

Bonacic-Koutecky et oiP 

MRD-CI 

0.63 

H-F 

-60.135(1) 

OweiP 1 

DMC 

0.67(2) 

PW91-GGA PW91-GGA -60.136(1) 

Nissenbaum et 1 

DMC 

0.733(4) 

PBE-GGA 

PBE-GGA -60.144(1) 

Rao et alW^ 

CI-SD 

0.7375 

PBE-GGA 

PZ-LDA -60.151(1) 

Wheeler et aZP 

CCSD(T) 

0.7445 



present work 

DMC 

0.744(3) 




expt. 

0.84(5) 



Brechignac et 

expt. 

0.63(4) 

ing hypothetical suggestion has been raised by Neaton 


involve states which would lie very high in energy since 
such excitations would require strongly localized type of 
orbitals. The fact that the nodal surfaces are minimally 
affected by the Is sub-shell is also suggested by the highly 
accurate calculations presented above for Li dimer. 

The table includes also experimental data from the two 
available sources, WtP3 and Brechignac et a/Fp 0.84(5) 
and 0.63(4) eV respectively. These data show signifi¬ 
cant differences and seem inconsistent with each other. 
Considering a reasonable agreement between the four in¬ 
dependent calculations in Table m we essentially claim 
that the experimental data is very inaccurate and our 
present calculation, being produced by an upper bound 
method and having the lowest total energy of the the¬ 
oretical calculations, is the most accurate estimation of 
this total energy to date. 


III.4. Bulk Lithium in a BCC crystal 

Because of its position as the bulk crystal with only 
one valence electron per atom, lithium in the body- 
centered-cubic (BCC) crystal (Pearson symbol cI2) has 
been stu died a few times by QMC methods in the liter- 
atu nj 15 b 6 * 43 l 5 H Surprisingly though, none of these cal¬ 
culations have used FN-DMC with the core electrons in¬ 
cluded. 

Further interest in studying lithium crystals with QMC 
methods was stimulated by recent experimental and the¬ 
oretical developments. For example, interesting phenom¬ 
ena for high pressures, including superconductivity, have 
been reported by experimental studies l 14 l 22 l . An intrigu- 


and Ashcroft that lithium solid may undergo a Peierls 
transition into a so-called “alkali electride” at high pres¬ 
sures, an exotic phase where paired Li ato ms are stabi¬ 
lized by pockets of highly localized electrons 31 ^ 

We treat lithium solid in the cI 2 crys tal by FN-DMC at 
the experimental lattice constant! 16 * 45 ! ao = 3.482 A. Our 
FN-DMC calculations use a single determinant Slater- 
Jastrow wave function with orbitals taken from DFT cal¬ 
culations. Since the nodal surface would be rather diffi¬ 
cult to improve upon within the QMC calculation of crys¬ 
talline system, we begin by comparing several DFT func¬ 
tionals to find the best nodal surface; and we found the 
results to lie within 2 — 2.5 mHa per atom of the highest 
quality nodal surface (PBE-PZ functional). This strat¬ 
egy has been motivated by our work on DFT gene rated 
orbital sets to find the most optimal nodal surfacd^^l. 
A select subset of the functionals tested are reported in 
Table VII The total energy is integrated over the irre¬ 
ducible Brillouin zone by the so-called “twist-averaging” 
procedure: a DMC calculation is carried out for each 


symmetry unique fc-point in a uniform 8x8x8 Monkhorst- 
Pack mesh and the resulting energies are weight-averaged 
using bv the geometric multiplicity of the k-point as the 
weighfpH To treat the finite size errors that occur, both 
due to the Ewald sums and the finite number of twists, 
we collect statistics on the static structure factor S(k) 
during the QMC simulations. These data are plotted in 
Figure [3] The correction to the finite-size errors in the 
simulation cell’s energy is calculated using the functional 
form for S(k ) as detailed by Chiesa et aZF 1 ®. The 8-atom 
(triangle symbols) and 16-atom (circle symbols) data for 
variational Monte Carlo (VMC) and DMC mixed esti¬ 
mators are less well converged when compared to the 54- 
atom cell. The corrected DMC mixed estimator (blue) 
is, however, consistent for all sizes of cell. This suggests 
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Simulation Cell Size (N 1 ) 


FIG. 3. QMC results for S(k) for several sizes of simula¬ 
tion cell. Square symbols denote results from the 54-atom 
simulation; circle symbols, 16-atom simulation; and triangle 
symbols, 8-atom simulation. The curves shown are fit to the 
54-atom data with the function S(k ) = exp{ — ak b }. For 
small k, RPA predicts 6 = 2. For the fits shown, the 6 param¬ 
eter for VMC are 1.45; DMC with the mixed estimator, 1.67; 
and the DMC extrapolated estimator, 1.99. 


that, at least for some systems, it is possible to estimate 
the static structure factor correction accurately with data 
from smaller simulation cells. Since the systems is a sim¬ 
ple metal, within the Random Phase Approximation the 
behavior of S(k) for small values of k is expected to be 
proportional to k 2 , as detail by Bohm and Pineal The 
curve fit to our calculated values indicates that S{k) in 
our simulation is oc k 199 so that we have reasonable con¬ 
fidence in the quality of the corrected mixed DMC esti¬ 
mator resultP20. We use a linear fit to the equation 

E n =E oa - ^ (6) 

to extrapolate the total energy per atom to the infinite 
bulk of the cI2 lithium crystal, as shown in Figure [4] 
The data in Figure [4] show that after twist averaging 
and corrections, the 16 atom cell energy per atom is 
ss 2.5 mHa from the infinite bulk value, and the 24- 
atom cell, ss 1 mHa, while the 54 atom is less than a 
mHa from the extrapolated value. We believe that the 
Li bulk energy could be further improved by more so¬ 
phisticated orbitals and optimization by employing more 
accurate wave functions, such as pfaffians 4 5 ; however, in 
this study the focus was to understand the trial functions 
that are, at present, commonly used for solid state and 
quantum chemistry calculations. 

III.5. Summary 

In order to compare the quality of results for these sys¬ 
tems, we estimate the exact energy per atom as the en¬ 
ergy of the single atom and the zero point motion (ZPM, 
infinite nuclear mass) corrected binding energy per atom, 


FIG. 4. The total energy per atom for twist averaged FN- 
DMC with and without the finite size error (FSE) correc- 
I iomP'* plotted against the inverse of the number of atoms in 
the cell. The data are fit with an extrapolation to infinite bulk 
size. The infinite bulk total energy per atom of -7.5371 Ha 
for the FN-DMC, and -7.5369 for FN-DMC with FSE correc¬ 
tions. These values are within error bars of each other. The 
statistical error bars on the data are smaller than the size of 
the plot symbol and so are not visible. 


TABLE VIII. Summary of the estimated total energies per 
atom for a sequence of different size Li systems. E to t for 
n — 4 cluster we substitute use the obtained FN-DMC value 
for the binding energy. 


size 

Efn-dmc 

Est. Eo Estimated from 

i 

-7.4780(1) 

—7.47806 Hylleraas expan. 

2 

-7.4976(lf 

1 -7.4977 Exp.+ZPlvfj 

4 

-7.50538(1] 

—7.50541 FN-DMC binding 


cI2 crystal -7.5369(6) -7.54078 Exp.+ZPk|] 

a from referencfUki 

b Atom + ZPM corrected exp. binding from Ref.]] 
c Atom + ZPM corrected exp. binding from Ref ! 20 ! 42 ] 


in a manner similar to the previous studies, see for ex¬ 
ample, RefP. For the Li 4 cluster, we substitute our own 
FN-DMC calculated binding energy as the best currently 
available estimation. Table |VTH| shows the per atom en¬ 
ergy evolving toward the bulk value as the size of the 
lithium system increases. 

The imperfect result for the FN-DMC correlation en¬ 
ergy in the Li solid is still very accurate, higher only by 
about 3 mHa (~ 3 %) from the estimated exact value 
obtained from the experimental cohesive energy. The 
underestimation of the cohesive energy is essentially the 
same magnitude, approximately 0.1 eV. We consider that 
remarkably accurate in light of how simple the single¬ 
reference Hartree-Fock wave function nodes are. Given 
the fact that the spin-up and spin-down channels are 
completely decoupled in HF, the complexity of the nodal 
surface is not fully captured by this trial function; nev- 
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ertheless the accuracy of the total energy appears to be 
quite robust so that the impact of these errors is compa¬ 
rably small. We conjecture that the electronic structure 
is dominated by the nearly-free electron picture that is 
far from the strongly correlated regime so that single¬ 
reference wave functions lead to qualitatively and also 
quantitatively accurate descriptions of ground state prop¬ 
erties of the Li solid. 


IV. CONCLUSION 

As the size of Li systems increases from a single atom 
to the bulk crystal, it is clear that the complexity of 
the nodal hypersurface grows. In the simplest case of 
the atom, a nearly exact approximation to the node is 
known. For related small systems, the nodal error er¬ 
rors are small, and it is possible to recover almost exact 
nodes with acceptable sizes of expansions in excited de¬ 
terminants. What is valuable and somewhat unexpected, 
is the fact that the accuracy of the FN-DMC calculation 


with single reference trial functions is high even for the 
Li solid. Note that the solid phase is metallic, so that its 
electronic structure is different from atomic and molecu¬ 
lar systems with localized ground states. The presented 
calculations show that for Li systems, readily available 
trial wave functions are sufficiently accurate to provide 
cohesive and binding energies to within an accuracy of 
0.05 — 0.1 eV. While our understanding of the fixed-node 
errors is still limited, the results presented here uncover 
another piece in the mosaic of previously obtained results 
which indicate that both the electronic density and the 
complexity of bonds, in particular, the bond multiplici¬ 
ties, strongly influence the nodal accuracy. 
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